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ABSTRACT 

We investigate the formation and evolution of a first core, protostar, and circumstellar 
disc with a three-dimensional non-ideal (including both Ohmic and ambipolar diffu¬ 
sion) radiation magnetohydrodynamics simulation. We found that the magnetic flux 
is largely removed by magnetic diffusion in the first core phase and that the plasma j3 
of the centre of the first core becomes large, /3 > 10'^. Thus, proper treatment of first 
core phase is crucial in investigating the formation of protostar and disc. On the other 
hand, in an ideal simulation, /3 ^ 10 at the centre of the first core. The simulations 
with magnetic diffusion show that the circumstellar disc forms at almost the same time 
of protostar formation even with a relatively strong initial magnetic field (the value for 
the initial mass-to-flux ratio of the cloud core relative to the critical value is /i = 4). 
The disc has a radius of r ~ 1 AU at the protostar formation epoch. We confirm 
that the disc is rotationally supported. We also show that the disc is massive {Q ^ 1) 
and that gravitational instability may play an important role in the subsequent disc 
evolution. 

Key words: star formation - circumstellar disc ~ methods: hydrodynamics - 
smoothed particle hydrodynamics - protoplanetary disc - planet formation 


1 INTRODUCTION 

The molecular cloud core is the f ormation site o f the star. Al¬ 
ready almost half-a-century ago, iLarso nl ll 19691) investigated 
the formation process of the protostar with one-dimensional 
radiation hydrodynamics simulation starting from a grav¬ 
itationally unstable cloud core. An overview of the evolu¬ 
tion obtained from that simulation is as follows: While the 
dust thermal emission effectively removes the thermal en¬ 
ergy generated by the compressional heating caused by the 
gravitational collapse, the gas evolves almost isothermally. 
At p ~ 10“^®g cm~®, the compressional heating overtakes 
the radiative cooling and the gas begins to evolve adiabat- 
ically. In this adiabatic evolution phase, the temperature 
evolves as T oc where 7 is the adiabatic index (7 = 5/3 
for r < 100 K and 7 = 7/5 for 100 < T < 2000 K). Be¬ 
cause this index is larger than the critical adiabatic index 
for spherical gravitational collapse, 7crit = 4/3, the gravi¬ 
tational collapse temporarily halts and a quasi-hydrostatic 
core forms, commonly known as the first core. When the 
central temperature of the first core reaches ~ 2000 K, the 
hydrogen molecules begin to dissociate. This endothermic 


reaction changes the effective adiabatic index to 7eff = 1.1. 
Because this is smaller than 7crit, the gravitational collapse 
resumes, which is known as the second collapse. Finally, 
when the molecular hydrogen is completely dissociated, the 
gas evolves adiabatically again and the gravitational collapse 
finishes. The adiabatic core formed at the centre is the pro¬ 
tostar (or the second core). This evolution process was later 
confirmed, more sophisticated one-dimensional simula tions 
llMasunaga fc Inutsukal[200^ : IVavtet et al.ll20l3 . l2013l) . 

Although, the general picture of the formation pro¬ 
cess of the protostar was established bv iLarso 3 (Il969l) with 
one-dimensional simulations, multidimensional simulations 
are necessary to investigate important phenomena such as 
the formation and evolution of the circumstellar dis c. After 
the ra diation hydrodynamics simulations done by iLarsonl 
lll969l) . it took several decades to develop and perform 
three-dimensional rad iation hydrodynamics simulations of 
gravi t ational collapse Iwhitehouse &: Batell2006l: lBatell2()10l. 


20U; Tomrdn^jt_aJj^i Tsukamoto. Machida fc Inutsukal 


201s ; Tsukamoto et aLll2015l) . These studies revealed that 
the multi-dimensionali ty ca u ses n ew and interesting phe¬ 
nomena. For example, iBatd ll2010l) found that the bipolar 




































2 Tsukamoto et al 


outflow from the first cor e can be driven by radiat ive feed¬ 
back from the protostar. iTsukamoto et aki (120151') investi¬ 
gated the evolution of the circumstellar discs in the unmag¬ 
netized cloud core and found that the temperature structure 
of the disc is determined by diffusive radiative transfer in the 
radial direction in its early evolution phase. 

The magnetic field is another important ingredi¬ 
ent in the star formation process. Observations sug¬ 
gest that the molecular c l oud cores are magnetize d 
fe.g. iHe iles fc Troland l2(105l : iTroland fc Crutchen I2OO8I I. 
IXroland &: Crutcher! ( 2008l l showed that the mean value 
of the mass-to-flux ratio relative to the critical value, /r, 
of the nearby dark cloud cores is /r ~ 2 — 3 and sug¬ 
gested that the magnetic field of the typical cloud core 
is relatively strong. The magnetic field drives the out¬ 
flow from both the first core and the protostar. The 
typical velocity of the outflow is determined by the ro¬ 
tational velocity at the launching point (u ~ 2 km/s 
from th e first core an d u ^ 20 km/s from the pro¬ 
tosta r) llToniisalM 200^ Machid^ Inutsuka fc Mat sumotc 


200s; iHennebelle fc Fromand bOOsi r iPrice. Tricco fc Bate 


2012f ). Another important effect caused by the magnetic field 


is the removal of the gas angu lar momentum. This effect 
is kn own as magnetic braking (iMouschovias fc Paleologoul 
[T^. Until recently, it was believed that the disc forma¬ 
tion is a natural consequence of the gravitational collapse of 
a rotating molecular cloud core. Actually, three-dimensional 
simulations, with a weak magnetic field or without it, show 
that a relatively large circumstellar disc (with a radius of 
several ten s of AU) de v elops i n the early phase of protostar 
formation ( 20n Tsulamoto^fc_Machidal l201ll . 


l2013l : [Tsukamoto. Machida fc Inutsukall2013l ') However, pre¬ 
vious works with ideal magnetohydrodynamics (MHD) sim¬ 
ulations have shown that the relatively strong magnetic field 
{fj, ~ 1) completely suppresses the formation of a rotation- 
ally su pported disc around the protostar at its forma tion 
epoch iMellon &: Lill200^ : IHennebelle fc FromandbOOSl l. 


Ideal MHD is, however, not a good approximation for 
the simulations of the magnetized molecular cloud core. Be¬ 
cause the ionization degree of the cloud core is quite low, it 
is expected that non-ideal magnetic effects such as Ohmic 
diffusion. Hall effect, and ambipolar diffusion play important 
roles during the formation and evolution of the circumstellar 
disc. 


The influence of non-ideal magnetic effects on the disc 
forma tion is still controversial. iLi, Krasnopolskv fc Shan3 
(l201lli investigated the influences of the non-ideal magnetic 
effects. They pointed out that ambipolar diffusion is the 
dominant diffusion process of the magnetic field and con¬ 
cluded that neither Ohmic nor ambipolar diffusion weak¬ 
ens the magnetic braking and that the disc formation is 
still strongly supp ressed even with the magnetic diffusion . 
On the other hand. lMachida. Inutsuka fc Matsumotol il201ll f 
showed that a relatively large disc of about a few tens of AU 
in size forms in the early phase of the protostar formation 
although they considered only Ohmic diffusion. 

The discrepancy could come from the difference in the 
initial conditions and the treatment of the inner boundary 
(or a si nk at the cen t re) of the simulations. In the simula- 
tions of iMellon fc Lil (l2008ll and lLi, Krasnopolskv fc Shangj 
(l201iri . the inner boundary or sink is set from the be¬ 
ginning of the simulations. In such a set-up, the simu¬ 


lations cannot follow the evolution of a first core which 
is mainly supported by gas pressure and not necessarily 
by rotation. Although the first core is a transient object, 
its density is high enough that the magnetic flux is ef- 
ficiently removed from the first core during its evolution 
dPapD. Basu fc Kunj|20l3 '). Furthermore, it is suggested 
that the greater p art of the first core directly b ecomes the 
circumstellar disc (iMachida fc Matsumotol 1201 il l just after 
the protostar formation. Therefore, calculating the first-core 
phase correctly in the simulations is crucial to investigate 
the very early phase of disc evolution. On the other hand, 
IMachida. Inutsuka fc Matsumotol 1I2OIIII used sink cells with 
“threshold density”. In their simulations, the sink cell takes 
in the gas when its density becomes larger than the threshold 
density. In this case, the gas whose density is smaller than 
the threshold density can reside inside or around the sink 
cell regardless of whether the gas is rotationally supported 
or no t. This treatment may also affect the disc evolution pro¬ 
cess. IMachida, Inutsuka fc Matsumot3 (I20l4) showed that 
the sink treatment (its radius and the threshold density) 
significantly affects the formation and evolution of the cir¬ 
cumstellar disc. 

To reveal the realistic formation and evolution processes 
of the first core, the protostar, and the circumstellar disc, ap¬ 
propriate treatment of the radiation transfer in the simula¬ 
tion is crucial, because the magnetic diffusion coefficients are 
functions of temperature. The previous studies with MHD 
simulations mentioned above do not include radiation trans¬ 
fer and employ a simplified equation of state (EOS) which 
mimics the temperature evolution of the centre of the cloud 
core. We call this the barotropic approximation. The simula¬ 
tions with radiation transfer, however, have shown that the 
temperature structures in the first core or around the pro¬ 
tostar are strikingly differe nt from those expected from the 
barotropic approximation (|Whitehouse_J^^atj 20061 : iBatd 
I2OI0I : iTomida et akll^lSl : ITsukamoto et al.ll201f ). 

Three-dimensional simulations which include both the 
magnetic field an d radiation transfer h ave not been success¬ 
ful until recently. iTomida et al.l (l2013l l was the first to suc¬ 
ceed with such a simulation with a grid code and found 
that the Ohmic diffusion alters the structure around the 
protostar significantly. With ideal radiation magnetohydro¬ 
dynamics (RMHD) simulations using the smoothed particle 
hydrodynamics (SPH) method, iBate. Tricco fc Prica (l2014l l 
also investigated the formation and evolution of the pro¬ 
tostar, especially the long-term evolution of the bipolar jets 
driven around the protostar. They showed that the jets heat 
up the gas in the envelope after they break up the remnant 
of the first core. Such a radiative heating process may affect 
the ionization degree of the g as and change the m a gnetic 
diffusion coefficients. However, iBate. Tricco fc Pried (l2014l l 
did not consider magnetic diffusion processes. 


As pointed out 


previous 


studies 


dLi, Krasnopolskv fc Shaji3 1201 il l, it is expected that 
ambipolar diffusion will play a role during the formation 
process of the protostar and the dis c arou nd it. Very 
recently, ITomida. Okuzumi fc Machidal (l2015ll conducted 
a simulation with both Ohmic and ambipolar diffusion. 
However, they only calculated the evolution until the end of 
the first core phase with ambipolar diffusion and the effect 
of the ambipolar diffusion is still unclear. 

In this paper, we investigate the formation of the first 
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core, protostar, and the circumstellar disc using a three- 
dimensional non-ideal RMHD simulation. We employ the 
SPH method and use it to produce the first results of the 
three-dimensional non-ideal RMHD simulations with SPH. 
Here, we focus on the effects of magnetic (Ohmic and am¬ 
bipolar) diffusion, but do not include the Hall effect. To 
avoid the numerical artefact caused by the sink, we do not 
introduce it, but rather investigate the structure around the 
protostar to determine whether the formation of the cir¬ 
cumstellar disc is possible at the very early phase of proto¬ 
star formation. This paper is organized as follows: In §2, we 
briefly describe the non-ideal magnetohydrodynamic effects. 
In §3, we describe the numerical method and initial condi¬ 
tions for the simulations, the results of which are given in 
§4, and then summarized and discussed in §5. 


2 NON-IDEAL MAGNETOHYDRODYNAMIC 
EFFECTS 


The ionization degree in the molecular cloud core is quite 
low and the gas can be regarded as weakly ionized plasma. 
In weakly ionized plasma, the microscopic collisions between 
neutral, positively-charged, and negatively-charged particles 
produce finite conductivity and non-ideal magnetohydrody¬ 
namic effects, or in short, non-ideal effects arise. 

The non-ideal effects appear as the correction terms in 
the induction equation. They can be derived by calculating 
the drift velocity of the charged particles. Here, we derive the 
in duction equati o n for the weakly ionized plasma according 
to lWardle fc Nd (1 19991 '! and I WardI3 (l2007l ~l. 

We start with 

= -cV X E, (1) 


dt 


J = —V X B. 

47r 


( 2 ) 


where B is the magnetic field, J is the current density, E is 
the electric field, and c is the speed of light. By the Lorentz 
transformation to the rest frame of the fluid (that is essen¬ 
tially the rest frame of bulk of neutral particles), the electric 
field becomes 

E' = E+^. (3) 


Here, v and E' are the fluid velocity and the electric field 
in the rest frame of the fluid, respectively. The conductivity 
in the weakly ionized plasma can be calculated using the 
balance of the force that acts on the charged particles. 


Zje(E' -b ^ - ^jprrijVj = 0. 


(4) 


Here, subscript j denotes the species of charged particles, 
ZjC is the charge, vj is the relative velocity of charged par¬ 
ticles in the fluid rest frame, 'yj = {av)j/{mj + m) where 
{(jv)j is the rate coefficient for momentum transfer, rrij is 
the mass of charged particles, m is the mean mass of neu¬ 
tral particles, and p is the density of neutral particles. Note 
that, in the weakly ionized plasma, most of the particles are 
neutral and the inertia of the charged particles and the col¬ 
lisions with other charged particles are negligible. Note also 
that, under the MHD approximation, the difference between 
the magnetic field and the current density in computation 
frame and those in the rest frame is negligible. We assumed 


the local charge neutrality "YijnjZj = 0. By inverting (|4} 


for Vj and calculating the current density, J = 
we obtain 

UjZjevj, 

J = (JoE^ -b (JjjB X E^ — (ap — (to)B X B X E\ 

(5) 

where 


= -g nj Zj Pj , 

(6) 

ec rij Zj 

" B^l + py 

3 

(7) 

ec UjZjPj 

p 1 + py 

3 

(8) 


(9) 


are the Ohmic, Hall, and Pedersen conductivities, respec¬ 
tively. Here, Pj = ZyeBjirrijC^jp') is the Hall parameter 
which is the product of the cyclotron frequency and stop¬ 
ping time. Finally, by inverting equation ([5| for E' and using 
equation o and m, we obtain 

f = 

— V X |»?o(V X B) -b riH{V X B) X B — rjA{{V x B) x B) 

This is the induction equation with non-ideal effects. The 
second, third, and fourth term on the right hand side of 
equation GQl describe the Ohmic diffusion. Hall term, and 
ambipolar diffusion, respectively. Here, 


VO — 

47r(Jo ’ 

Vh = 

an 

+ af) 

VA = 

c?<jp 

47 r(cr| -b cr|,) 


( 11 ) 

( 12 ) 

(13) 


are the Ohmic, Hall, and ambipolar diffusion coefficients, 
respectively. In this paper, the Hall term is neglected owing 
to the numerical difficulty associated with it. The effect of 
the Hall term will be investigated in future works. 

We constructed the data table of the diffusion co¬ 
efficients by calculating a chemical reaction network of 
H^, HCO*^, Mg'*', He"*", C"*", H'*', e“ in gas phase and the 
positively-charged, neutral, and negatively-charged dust 
grain of uniform size using the meth ods descr i bed i n 
iNakano. Nishi fc Umebavashif ll2002h and lOkuzurnU (l2009ll . 
We assumed that the dust to gas ratio is 10“^. We also 
assumed that the dust grain size and density are a = 
3.5 X 10“^ pm and pd = 2 g cm“®, respectively. We con¬ 
sidered non-thermal ionization by the cosmic rays and ther¬ 
mal ionization in our calculations. The cosmic-ray ionization 
rate was fixed to be ^cr = When the tempera¬ 

ture reaches T ~ 1000 K, thermal ionization is the dominant 
source of ionization. In this paper, we consider the effect of 
the thermal ionization by considering the thermal ionization 
of potassium. The coupling between the magnetic field and 
the gas quickly recovers around T ~ 1000 K because the 
thermal ionization provides a sufficient ionization degree. 

In figure [T] we show the Ohmic and ambipolar diffusion 
coefficients under the typical evolution of the gas. To make 
figure [T] we assumed that the temperature and magnetic 
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field change as, 

/ s 2/3 

The figure shows that the diffusion coefficients suddenly 
drop around p = 5 x 10“® g cm~® where the temperature 
is about T = 1000 K and the ionization degree quickly in¬ 
creases owing to the thermal ionization of potassium. 


3 NUMERICAL METHOD AND INITIAL 
CONDITIONS 


In this study, we solve the non-ideal radiation magnetohy¬ 
drodynamics equations with self-gravity. 


Dt 


D 

Dt 

_D 

Wt 


Dv 

1 /' 
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V p 
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-V X 


Er 
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V^<I> 
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• V V 


(16) 


V-F 


1 , 


|»7o(V X B) — p^(V X B) X B X b| , 
-b Kpc{arTg - Er), (17) 


1 


= —V-(P-b-B^)v-B(B-V 

P I 2 

— Kpc(arTg — Er) — V • Vd* 

— iv-[{(r/o(V X B) 

P 

— ?7a(V X B) X B X B 
= 'EvGp. 


}xb]. 


(18) 

(19) 


Here, p is the gas density, v is the velocity, B is the magnetic 
field, B is the unit directional vector of the magnetic field, 
P is the gas pressure, Er is the radiation energy, Fr is the 
radiation flux, is the radiation pressure, Tg is the gas 
temperature, Kp is the Plank mean opacity, e = pu+i(pv^-|- 
B= ) is the total energy with u specific internal energy, and 
$ is the gravitational potential. Parameters, Or and G are 
the radiation and gravitational constants, respectively. 

We adopt the gray approximation (frequency-integrated 
radiation transfer) and we assume local thermodynamic 
equilibrium (LTE). To close the equations for radiation 
transfer, we employ flux-limited diffusion (FLD) approxi¬ 
mations. 


Fr 

R 


n 




IVEr 


KRpEr ’ 

2 2 

VEr 


= OEr, 


n (g) n, X = ^ + ^ R^, 


|VEr| ■ 


Here, kr is the Rosseland mean opacity. FLD is a dif¬ 
fusion scheme which is designed to maintain the causal¬ 
ity of |Fr| < cEr. It is suitable for optically thick gas 


owing to its diffusive nature. In this paper, we use the 
SPH method to investigate the formation of a protostar 
and disc. The SPH method can be easily implemented 
and is suitable for simulations which treat the large dy¬ 
namic range because of its adaptive nature. The ideal MHD 
part was solved by adopting the Godunov smoothed par¬ 
ticle magnetohydrodynamics (GSPMHD) method in which 
the Godunov method and the method of characteristics 
are used to calculate the interactions b etween the particles 
inste ad of artificial dissipation terms (llwasaki fc Inutsukal 
I 2 OIIII . The divergence-free constraint on the magnetic field 
was maintained using the hyperbolic divergence cleaning 
method for GSPMHD (llwasaki fc Inutsukal 1201^ . The ra- 
diative transfer was treated by the FLD-SPH method 
jV^itehouse fc Batelf2003 : IWhitehouse, Bate fc MonaghanI 
I 2 OO 5 II . We treated Ohmic and ambipolar diffusion with 
the m ethod described by Tsukamoto. Iwasaki fc Inutsukal 
1 I 2 OI 3 II and IWurster. Price fc Avliffel ( 2014l ~). respectively. 
Both diffusion processes were accelerated by super t ime step¬ 
ping method (lAlexiades. Amiez fc Gremaudlll996h . To cal¬ 
culate the self-gravity, we adopted the Barnes-Hut tree algo - 
rithm with opening angle of 0 — 0.5 dBarnes fc Hud[i986l l. 
We do not use the individual time-steps and the particles 
are updated simultaneously. 

We ado p ted the tabulated EOS used in 
iTomida et al.l (l2013l ~). in which the internal degrees 
of freedom and chemical reactions of seven species 
Ha, H, H+, He, He+,He++,e" are included. We assumed 
that the hydrogen and helium mass fractions were X = 0.7 
and Y = 0. 28, respectively. The d ust opacity table was ob¬ 
tained from ISemenov et al.l ll20 03ll and the g as opacity table 
was obtained from Ferguson et al.l (l2005h . The resistive 
model is described in §2. 

We modelled the initial cloud core using an isothermal 
uniform gas sphere. The initial mass and temperature of 
the core were fixed at 1 M© and 10 K, respectively, with 
an initial core radius of i? ~ 3.0 x 10® AU. The core is ini¬ 
tially rigidly rotating with an angular velocity of Dq = 2.2 x 
10~®® s“®. The product of the angular velocity and the free- 
fall time is tfiDo = 0.19. The initial magnetic field was uni¬ 
form and parallel to the rotation {z-) axis with a strength of 
Bo = 1-7X 10^/tG. The corresponding initial mass-to-flux ra¬ 
tio relative to the critical value was p, = {M/^)/{M/^)cTit = 
4 where d* = ttR^Bq. We adopted a critical mass-to- 
fiux ratio of (M/<I>)crit = (0 .53/37r)(5/G)®^^ suggested by 
iMouschovias fc Snitzerl (Il976l l. The initial cores were mod¬ 
elled with about 10^ SPH particles. The boundary condi¬ 
tions of radiation transfer were introduced by fixing the gas 
temperature to be 10 K when the gas density was less than 
2.0 X 10“^^ g cm"®. 

We performed three simulations with and without 
Ohmic and ambipolar diffusion. The model names and the 
diffusion processes included in each model are summarized 
in Table 1. 


4 SIMULATION RESULTS 

4.1 Evolution at the centre of the cloud core 

To investigate how the magnetic field evolves during the 
gravitational collapse, we show the evolution of the central 
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Table 1. The model names and the magnetic diffusion they include. “Yes” means that the corresponding magnetic diffusion is included 
while a “No” means that it is not. 


Model 

Ohmic diffusion 

Ambipolar diffusion 

1 

No 

No 

2 

Yes 

No 

3 

Yes 

Yes 



Figure 1. Diffusion coefficients, rjo (solid) and rjA (dotted) as 
a function of density. For this plot, we assumed that the temper¬ 
ature and magnetic field are functions of density (see eq. Gl). 


magnetic field as a function of the central density in figure 
[21 At first, the magnetic field evolves as Be oc pc^^- This 
evolution is expected from a spherically symmetric collapse 
during which the central magnetic field evolves as Be oc R~^ 
due to the conservation of the magnetic flux, where R is the 
radius of the cloud. On the other hand, the central den¬ 
sity evolves as pe oc R~^ or, equivalently, R oc pO^^ ■ Thus, 
Be oc R~^ oc pc^^■ The increase in the magnetic field almost 
stops {Be oc Pe) at 10“^® ^ Pc ^ 10“^^ g cm“® because the 
Lorentz force becomes comparable to the gravitational force 
and the gas moves almost parallel to the magnetic field. The 
a-component of the magnetic field still dominates other com¬ 
ponents and the gas moves almost vertically. As a result, a 
sheet-like structure (or pseudo disc) forms. When the cen¬ 
tral density reaches pe ~ 10“^® g cm~®, the central magnetic 
field evolves as Be oc pp which indicates that the collapse 
becomes sheet-like. In the gravitationally collapsing isother¬ 
mal sheet (whose scale-height \s H — cf/ (GE) = CsjOGpc ), 
the central magnetic field and density evolves as Be oc R~^ 
and Pe oc R~^H~^ oc R~^ and hence Be oc pV^■ 

Once the central density exceeds p ~ 10“^^ g cm~®, 
the magnetic diffusions becomes effective and the magnetic 
freezing is no longer valid for resistive models. The mag¬ 
netic flux is removed from the central part and the difference 
between the ideal model and resistive models can be seen. 
The central magnetic field of model 1 (the ideal model) is 
about sixty times larger than that of model 3 (with Ohmic 
and ambipolar diffusion) when the central density reaches 
Pc ~ 10“® g cm“®. Around the pe ~ 10~® g cm“®, the mag¬ 
netic diffusion becomes ineffective again owing to the ther¬ 



Figure 2. The evolution of the central magnetic field as a func¬ 
tion of central density. The solid, dashed, and dotted lines show 
the results of model 1 (ideal), model 2 (with Ohmic diffusion), 
and model 3 (with Ohmic and ambipolar diffusion), respectively. 
Differences between the models can be seen when the central den¬ 
sity exceeds pe > g cm“® and magnetic diffusion becomes 

effective. 

mal ionization and the flux freezing recovers in the resistive 
models. This causes Be oc pV^ again. 

4.2 Structure of the first core 

When the central density reaches pe ~ 10“^® g cm“®, the 
gas becomes opaque and the compressional heating due to 
the gravitational contraction cannot radiate away efficiently. 
As a result, the gas evolves adiabatically and a pressure- 
supported core, the first core, forms. The first core phase 
lasts until the central temperature becomes Te ~ 2000 K 
(or pe ~ 10“® g cm~®) at which point the dissociation of 
hydrogen molecules begins. The durations of the first core 
phase are about 620 years for model 1, 810 years for model 
2, and 940 years for model 3. The first core phase is defined 
as the phase in which the central density is g cm“® < 

pe < 10“® g cm“®. The difference in the duration is due to 
the rotation of the first core. 

To investigate the structure in and around the first 
core, we show the cross sections of the density, gas tem¬ 
perature, and plasma j3 around the first core in the y = Q 
plane in figure [3] at the end of the first core phase {pe ~ 
3x 10“® g cm“®). The plasma /3 is defined as /3 = Pgas/TVnag 
where Pgas and Pmag are the gas pressure and magnetic pres¬ 
sure, respectively. Note that the box size of the density cross 
sections is four times that of the other cross sections to com¬ 
pare the outflow structures of each model. 

To obtain the cross section and the profiles, the phys- 
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x(AU) 




x(AU) 


x{AU) 


Figure 3. The cross sections of the density, gas temperature, and plasma (3 (from left to right) around the first core in the y = 0 plane. 
The top row corresponds to model 1, the middle row to model 2, and the bottom row to model 3. The thin black lines show the contour 
of each quantity, while the thick black lines in the density cross sections show the \vz \ = 0 contour. This traces the outflow regions. The 
red arrows in the density cross sections show the velocity field. The box size of the density cross sections is four times larger than the 
other cross sections to show the outflow structures. The color bars of the density, temperature, and plasma /3 show log(p [ g cm~^]), 
log(T [K]), and log(/3), respectively. 


ical quantities are needed at grid points. In this paper, the 
physical quantities are calculated at grid points through, 


/(Xgrid) 


'^j "T-J - Xj, hj ) 


( 20 ) 


In the left panels, we show the density cross section. The 
thick black solid lines show the Vz = 0 contour and trace 
the outflow structure. The outflow formed in both model 
1 and 2, but did not in model 3 at the epochs. Although 
the outflow did not form in model 3, we confirmed that the 
outflow does form in a simulation with both Ohmic and 
ambipolar diffusion when the initial rotation of the cloud 
core is slightly larger than in the model 3. Therefore, we 
conclude that the magnetic resistivity delays the formation 
of the outflow rather than suppressing it. In our results, 
both the magneto-centrifugal force and the magnetic pres¬ 
sure play a role in driving the outflow. 


In the middle panels, we show the temperature cross 
section around the first core. The high temperature (T ~ 
1000 K) regions with radius of r ~ 5 AU are formed at 
the centre due to the radiative transfer. The high tempera¬ 
ture region is extended compared to the case in which the 
barotropic EOS is adopted. Because the thermal ionization 
becomes effective at T ~ 1000 K, the coupling between the 


magnetic field and the gas recovers in the relatively large 
part of the first core when the radiative transfer is taken 
into account. This recoupling causes the amplification of the 
magnetic field inside the first core due to the rotation. 

In the right panels, we show the cross section of the 
plasma /3. Because of the magnetic diffusion, the magnetic 
flux is efficiently removed from the first core in the resistive 
models. Thus, in the resistive models, (3 > 10® at the centre 
of the first core while in the ideal model, jH ~ 10. After the 
removal of the magnetic flux, the coupling between the gas 
and the magnetic field recovers at the central region of the 
first core owing to the thermal ionization and the magnetic 
field in the first core is reamplified by the rotation. As a 
result, the plasma /3 around the centre slightly decreases in 
the resistive models. This amplification is clearly seen in the 
middle right panel. 

Figure 0] shows the profiles of the density, gas temper¬ 
ature, and plasma p at the same epoch of figure [S] In all 
models, the central density and the central temperature of 
the first core are pc ~ 3 x 10~® g cm“® and Tc ~ 10® K, re¬ 
spectively. The density and temperature profiles show that 
the first cores formed in each model have very similar struc¬ 
tures. This is because the angular momenta of the first cores 
are not significantly different and the structural difference 
caused by rotation is negligible. The density on the a:-axis 
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is larger than that on the 2 -axis outside of the first core be¬ 
cause the pseudo disc has formed in the x direction. On the 
other hand, the temperature profiles along the x and 2 -axis 
do not differ significantly and the temperature structure is 
hence almost spherically symmetric. 

Due to the magnetic diffusions, the plasma /3 in the 
central region of the first core differs significantly between 
the ideal model and the resistive models. In model 1, the 
plasma P inside the first core is ~ 10 and almost constant 
in the x direction. In the model 2, the plasma /3 at the centre 
of the first core becomes /i ~ 6 x 10®. This is hence about 
three orders magnitude greater than for the ideal model. The 
magnetic flux removed from the first core piles up around 
it and the plasma /3 on the x-axis becomes smaller than the 
ideal model at the perimeter of the first core [x ~ 10 AU). In 
model 3, the plasma P at the centre of the first core becomes 
/3 ~ 6 X 10“^, which is much higher than for the model 2. In 
the z direction, the plasma p quickly decreases in all models 
because of the large density gradient in this direction and 
the magnetic field amplification by the first core rotation. 
Because the plasma P is larger than 10 inside the first core, 
the magnetic pressure does not affect the pressure support 
in the first core. 

A notable difference between models 2 and 3 is the 
plasma P in the x direction at the perimeter of the first 
core. In the model 2, only Ohmic diffusion is considered. 
The Ohmic diffusion coefficient is an increasing function of 
density and does not depend on the magnetic field. Roughly 
speaking, the Oh mic diffusion does not play a role whe n 
p < 10“^® g cm“® (iMachida, Inutsuka fc Matsumot^l2007li . 
Because the density of the first core is p > 10“^® g cm“®, 
the magnetic flux piles up outside the first core. By this 
pile-up, the plasma P beyond the first core in model 2 is 
/? ~ 1 around x = 10 AU and becomes much smaller than 
for model 1 at larger x. In model 3, the ambipolar diffusion 
is included as well. The diffusion coefficient of the ambipolar 
diffusion is a function of the magnetic field and does not de¬ 
pend strongly on the density. Therefore, it is expected that 
the pile-up of the magnetic flux is smoothed by the ambipo¬ 
lar diffusion. Actually, the region of small plasma /3 (/I ~ 1) 
in the x direction broadens in the right panel. This difference 
can also be seen in the right panels of figured] 

In figure [S] we show the infall and rotation velocities 
along the x-axis and the infall velocity along the 2 -axis. The 
infall velocity along the a;-axis is larger than that along the 
2 -axis and the density on the i-axis is also much higher 
than on the 2 -axis at the surface of the first core (x, 2 ~ 10 
AU). Therefore, the mass accretion onto the first core is 
asymmetric and is maximal in the horizontal direction. 

The rotation velocity vp reaches its maximum value at 
a; ~ 2 AU in models 2 and 3. Inside this radius, the velocity 
profile obeys the rigid rotation relation, oc x. Note that 
a rigid rotation is expected when the density is constant 
because oc y' Mr/r oc por^ jr oc r, where Mr and po 
are the mass inside r and the density, respectively. In model 
2, Vcf, sharply decreases at the r ~ 4 AU. This is caused by 
the strong magnetic braking by the piled-up magnetic field. 
As mentioned above, the magnetic flux piles up around the 
first core. Hence, the magnetic braking is locally enhanced 
at r ~ 4 AU and the rotation velocity is decreases. The 
profile of the model 1 also obeys the relation of oc a; for 
X < 1 AU. On the other hand, for 3 AU < a; < 10 AU, the 


profile has a complex structure. This structure is also caused 
by the magnetic braking. Note that the plasma /3 inside the 
first core is still small in model 1. 

In figure [6| we show the evolution of the angular mo¬ 
mentum of the first core in relation to the central density. We 
define the first core as the region where p > 10“^® g cm~®. 
As we have seen above, the magnetic field in the first core 
becomes weak due to the magnetic diffusion which causes 
an inefficient angular momentum transfer by the magnetic 
braking. Thus, it is expected that the angular momentum 
of the first core becomes large in resistive models and it in¬ 
deed becomes large when the magnetic diffusion is included. 
The difference in the angular momentum between model 1 
and model 3 is a factor of 6 and hence, insignificant. Most 
of the initial angular momentum of the fluid element has 
already been removed during the isothermal collapse phase. 
With the parameters adopted in our simulations, a disc of 
r ~ 100 AU forms when the magn e tic field is neglected (see , 
e.g., iTsukamoto fc Machida! [io 111 : iTsukamoto et al.l 120151 1. 
Therefore, we conclude that the angular momentum of the 
first core depend more strongly o n the initial condition of the 
mole cular cloud cores (see, e.g., Ijoos. Hennebelle fc CiaiMil 

l2012tl . 

4.3 Formation of the protostar 

When the central density reaches pc ~ 10“® g cm“® and 
the hydrogen molecules are completely dissociated, the gas 
evolves adiabatically and the protostar forms at the centre 
of the first core. In figure [T] we show the cross sections of 
density, temperature, and plasma P around the protostar. 
The central density is pc ~ 10“® g cm“® at this epoch and 
just after the protostar formation. Note that the x, y, and 
color-bar scales differ between the ideal model and resistive 
models because the structure around the protostar in the 
ideal model is quantitatively different from the one in the 
other models. The density distributions of the resistive mod¬ 
els (middle and bottom left panels) exhibit the dumbbell-like 
structures. These structures indicate that the rotation plays 
a role in the resistive models. On the other hand, in model 1 
(the ideal model), the density structure is elliptical and there 
is no dumbbell-like structure even in vicinity of the proto¬ 
star. As we will show below, the rotationally supported disc 
quickly forms during the subsequent evolution in the resis¬ 
tive models but does not form in the ideal model. The tem¬ 
perature distributions around the protostar are smooth and 
roughly spherically symmetric in all models. The tempera¬ 
ture exceeds 1000 K and the magnetic diffusion is no longer 
effective in the entire region. In the model 2, the low P re¬ 
gion forms in the vertical direction. This structure is created 
by the rotational amplification of the magnetic field. As a 
result, the plasma /3 becomes P ~ 10“^. The magnetic field 
is also magnified in model 3. However, it is not a significant 
magnification and the plasma P in the vertical direction is 
still P ~ 10^ at this epoch. We cannot find any signature 
of the rotational amplification in model 1. The low P re¬ 
gion in the vertical direction is created by a dragging of the 
poloidal magnetic field. The figure [7] shows that the struc¬ 
tures around the protostar are significantly different even 
just after the protostar formation when the magnetic diffu¬ 
sion is considered. 

After the protostar forms, it evolves via the mass ac- 
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Figure 4. The density (left), gas temperature (middle), and plasma /3 (right) profiles. The epochs are the same as in the figure [3] The 
solid and dashed lines show the profiles of the x and z directions, respectively. The red, green, and blue lines show the results of model 
1 (ideal model), model 2 (with Ohmic diffusion), and model 3 (with Ohmic and ambipolar diffusion), respectively. 




Figure 5. The profiles of the infall velocity (left) and rotation velocity (middle) in the x direction and the infall velocity in the z 
direction (right). The epochs are the same as in the figure [3] The red, green, and blue lines show the results of model 1 (ideal model), 
model 2 (with Ohmic diffusion), and model 3 (with Ohmic and ambipolar diffusion), respectively. 


cretion from the remnant of the first core. In figure [SI we 
show the density and gas temperature along the x-axis 
(solid lines) and 2 -axis (dashed lines) at the end of the 
simulations. The central densities and temperatures reach 
Pc ~ 10“^ — 10“^ g cm“® and Tc > 10^ K, respectively. 
From the decrease in the density and temperature of the 
red lines around x,z ^ 10~^AU, we can identify the ra¬ 
dius of the protostar in the ideal model as r ~ 10“^AU. In 
the ideal model, the difference between the density in the 
horizontal and the vertical directions is not large and the 
density structure is almost spherically symmetric. On the 
other hand, the density profiles of the resistive models show 
a different structure around the protostar. After the forma¬ 
tion of the protostar, the rotationally supported disc of size 
1 AU quickly forms in resistive models in these epochs. Be¬ 
cause of the disc formation, the boundary of the protostar 
becomes ambiguous in the density and temperature profiles 
in the horizontal direction. Weak shock wave structures can 
be seen at a; ~ 1 AU in the green and blue solid lines of 
density. This is the boundary of the circumstellar discs. 

In figure [9l we show the infall and rotation velocity 
along the a;-axis. The left panel shows the infall velocity. In 
the ideal model, the infall reaches x ~ 10~^ AU, which shows 
that the first core remnant accretes directly onto the central 
protostar. On the other hand, the infall stops at a; ~ 1 AU in 
the resistive models. This radius corresponds to the shocks 
in the density profiles and thus to the edges of the discs. 
Note that there are the other shocks at a: ~ 10 AU. These 
are the accretion shocks at the surface of the first core. The 
remnant of the first core still exists in these epochs. 

We can see a clear transition of the rotation profile at 



Figure 6. The evolution of the angular momentum of the first 
core in relation to the central density. The solid, dashed, and 
dotted lines show the results of model 1, 2, and 3, respectively. 

X ~ 10”^ AU in the resistive models (blue and green lines). 
In a; < 10“^ AU, the profile obeys cc x and the gas rigidly 
rotates. This rigidly rotating region is the protostar and its 
radius in the resistive models is also r ~ 10~^ AU. In 10” < 
a; < 1 AU, the profile follows oc This is the rotation 

profile of the disc around the protostar. The rotation profile 
of the disc is more shallow than for a Keplerian disc (or disc 
subjected to a gravitational potential created by a point 
mass) which obeys the profile of Vr/, oc x~^'^. This means 
that both the self-gravity of the disc and the gravity of the 
central protostar influence the rotation profile. 
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4.4 Rotationally supported disc around protostar 


As we have seen above, there are several features of the 
density and velocity profiles which suggest the existence of 
a circumstellar disc. For example, the rotational velocity at 
the mid-plane of models 2 and 3 is considerably larger than 
the radial velocity in 10“^ ^ 2 ; < 1 AU. In addition, shocks 
exist at r ~ 1 AU in the density and infall velocity profiles. 
However, it is not clear from the above analysis whether the 
disc is rotationally supported or not. 

To confirm that the disc is really rotationally supported, 
the ratio of the sum of the centrifugal and the pressure gra¬ 
dient forces to the radial gravitational force, 


91 


vllr + VrPjp 
V,<E> 


( 21 ) 


is plotted in figure [TO] with the solid lines and the ratio of 
the centrifugal to the radial gravitational force, 


92 


\VlO\ 


( 22 ) 


with the dashed lines. Here, p and $ are the pressure and 
the gravitational potential, respectively. When gi = 1 and 
<12 qi, the gas is supported by the pressure gradient force. 
On the other hand, when gi = 1 and q 2 ~ 91 , the gas is 
mainly supported by the centrifugal force. 

The red lines show that ~ 1 and 52 ^ 9 i for 
X < 10~^AU. This means that a pressure supported sec¬ 
ond core (the protostar), whose radius is r ~ 10“^AU exists 
at the centre. On the other hand, the radial gravitational 
force always dominates other forces for 10~^AU < a; < 5AU. 
Therefore, neither the pressure gradient force nor the cen¬ 
trifugal force can cancel the gravitational collapse and no 
rotationally supported disc forms in the ideal model. On the 
other hand, the green and blue lines show that qi is almost 
unity for x < lAU and the gravitational force is cancelled 
in this region. For x < 10“^AU, the gi ~ 1 and 52 <C 91 , 
which shows the existence of a pressure supported protostar. 
Meanwhile, 52 is about 0.6 for 10“^AU ^ x < lAU and 60% 
of the gravitational force is cancelled by the centrifugal force 
and the remaining 40% is cancelled by the pressure gradient 
force in this region. Thus, the gas is supported mainly by 
the centrifugal force. From these results, we conclude that 
the rotationally supported disc forms naturally in the very 
early phase of the protostar formation when the magnetic 
resistivity is included and the first core phase is considered 
correctly. Note that the dips of the green and blue solid lines 
at the edge of the disc are due to the large pressure gradi¬ 
ent there. The ram pressure caused by the mass accretion 
should balance this. 

The first core directly becomes the disc and its mass 
is much larger than that of the protostar during its forma¬ 
tion epoch. Thus, it is expected that the self gravity plays 
an important role in the early phase of the disc evolution 
ilnutsuka, Machida fc Matsumotcj l2010l i. In figure 1111 we 
show Toomre’s Q value of the disc Q = 7 ^! where we 
approximate the epicycle frequency k as Q. In the disc re¬ 
gion 10“^AU < 2 ; < lAU Toomre’s Q value is Q ~ 2 — 3. 
As pointed out in previous studies, the disc becomes unsta¬ 
ble against non-axisymmetr ic perturbations when Q ~ 1.5 
and t he spiral arms develop (iLaughlin. Korchagin fc Adamj 
Il998h . The spiral arms invoke an angular momentum trans¬ 


fer. Although, the Q value is still slightly larger than 1.5, 
it is expected that the gravitational instability plays a very 
important role for the angular momentum transfer in the 
subsequent disc evolution because a large amount of the 
remnant of the first core is still accreating to the disc and 
the disc mass increases quickly. 


5 SUMMARY AND DISCUSSIONS 

In this paper, we investigated the formation and evolution 
of the first core, the protostar and the disc around the proto¬ 
star by using three-dimensional simulations with radiation 
transfer, as well as Ohmic and ambipolar diffusions. 

Our findings are summarized as follows. 

(i) The magnetic flux is largely removed in the first core 
phase. As a result, at the centre of the first core, plasma P 
becomes p > lO’^. On the other hand, the P at the centre of 
the first core in ideal simulation is /3 ~ 10 . 

(ii) Even though the plasma P inside the first core is 
significantly different in the resistive and the ideal mod¬ 
els, the angular momentum of the first core is not (within 
an order of magnitude). This is because most of the 
angular momentum has been removed before the mag- 
netic diffusion processes play a role. A ctually, figure 11 of 
Machida. Inutsuka fc Matsumot3 ll200it l suggests that most 
of the angular momentum is removed from the gas during 
the isothermal collapse phase. When the magnetic field is 
neglected, a disc with r ~ 100 AU forms in the cloud core 
for the parameters adopted in our simulations (se e, e.g., 
iTsukamoto fc Machidal 120111 : iTsukamoto et ahl l2015l l . This 
also suggests that most of the angular momentum is removed 
during the isothermal phase. 

(iii) With magnetic diffusions, a circumstellar disc forms 
around the protostar just after protostar formation even 
with a relatively strong initial magnetic field (we employ a 
uniform density sphere and an initial mass-to-fiux ratio rel¬ 
ative to the critical value of /r = 4). We confirmed that the 
disc is rotationally supported. The disc is massive enough to 
enable gravitational instability to develop in the subsequent 
disc evolution. Thus, the gravitational instability plays an 
important role in the early evolution of the circumstellar 
discs. 

The reason why most of the angular momentum is re¬ 
moved from the gas in the isothermal collapse phase can be 
understood by comparing the magnetic braking timescale 
fb ~ Aj/ua to the free-fall timescale tg, where Aj and va are 
the Jeans length and Alfven velocity, respectively. The mag¬ 
netic braking timescale is estimated as the time in which the 
inertia of the central region is equal to the inertia of the enve¬ 
lope where the Alfven wave sweeps dMatsumoto fc Tomisaka] 
l2004h . The ratio of the two timescale tb/tg is given as 
fb/fg ~ Aj/(uAtg) ~ y/P- In our simulations, the plasma P 
is p = 1.7 at the initial condition (p = 5.5 x 10“^® g cm“®) 
and decreases during the early isothermal collapse phase as 
P oc (?slv\ oc where we assume that Cs is constant 

and B oc p^^® as shown in figured When the central den¬ 
sity reaches pc = 10“^® g cm“®, tb/ig = \fP = 0.71 and the 
magnetic braking timescale becomes shorter than the free 
fall timescale. Therefore, tb/tg < 1 and the angular mo- 
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Figure 7. The cross sections of the density, gas temperature, and plasma 0 (from left to right) around the protostar in the J/ = 0 plane. 
The top, middle, and bottom row corresponds to model 1, 2, and 3, respectively. The thin black lines show the contours of each quantity. 
The color bars of the density, temperature, and plasma 0 are expressed as log(p [ g cm“®]), log(T [K]), and log(/3), respectively. Note 
that the x, y, and color-bar scales differ between the ideal model and the resistive models. 




Figure 8. The profiles of the density (left) and gas temperature (right) at the end of the simulations. The solid and dashed lines show 
the profiles in the x and 2 : directions, respectively. The red, green, and blue lines show the results of model 1, 2, and 3, respectively. 


mentum is largely removed during the isothermal collapse 
phase. 

Our results about the disc formation are largely consis¬ 
tent with those of the previous studies which followed the 
protostar formation with suf ficient resolution and consid¬ 
ered the first core p hase (e.g.. iMachida fc Matsumotoll20lll : 
IXomida et al.l 1201,31 ). We believe that the development of 
a disc at the very early phase of the star formation is a 
robust consequence. The previous research we mentioned 
above considered only Ohmic diffusion. On the other hand, 
we also included ambipolar diffusion. This does not change 


the overall formation process of the disc significantly. How¬ 
ever, it is possible that the ambipolar diffusion plays a more 
important role in the subsequent evolution of the disc be¬ 
cause it extends the density range in which the magnetic 
field and the gas are decoupled and allows the magnetic flux 
to escape from the disc. 

The difference in disc formation between the ideal 
model and resistive models is due to the strength of the 
magnetic field and not the difference in the angular momen¬ 
tum of the first core. In our simulations, the circumstellar 
disc forms in the resistive models (model 2 and 3) and does 
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Figure 9. The profiles of the infall velocity (left) and rotation velocity (right) 
[8] The red, green, and blue lines are also defined as in figure |8] 


in the x direction. The epochs are the same as in figure 


not in the ideal model (model 1). As we have seen above, in 
resistive models, the plasma /3 of the envelope around the 
protostar is /3 > 10^ except for the vicinity of the protostar 
of model 2 (the middle and bottom right panels of figure 
[7]l and the magnetic braking is ineffective. On the other 
hand, the magnetic field removes the angular momentum 
from the gas during the second collapse in the ideal model 
because the plasma /3 of envelope is 10“^ < /3 < 10^ (see, 
the top right panel of figure 0 and the magnetic braking 
timescale is comparable or less than the free-fall timescale 
(fb/iff ~ '/P)- This is why the circumstellar disc does not 
form in the ideal model. The simulation wit h Ohmic dif¬ 
fusion in iTomida. Okuzumi fc Machidal (|201!Th showed that 
the circumstellar disc forms even in the slowly rotating first 
core (J ~ 2 X 10®° g cm^ s“^ where J is the angular mo¬ 
mentum). Thus, the several-fold difference in the angular 
momentum does not affect whether or not the disc forms. 

Because the magnetic flux is largely removed in the 
first core phase, the proper treatment of the first core is 
necessary to investigate the formation of the protostar and 
disc. In previous works which argue that the disc forma- 
tion is strongly suppressed by the magnetic braki ng (e.g., 
iMellon fc Lil 20081 : iLi, Krasnopolskv fc Shanejl201ll ). the in¬ 
ner boundary was set from the beginning of the simulations. 
With this treatment, the previous works cannot follow the 
first core phase properly that should be supported by gas 
pressure. The discrepancy between our results aud those of 
these works should be du e to the different treatment s of the 
first core phase (see, also lPapp. Basu fc Kunj|2012l) . 

It is expected that the disc size becomes larger than the 
size obtained in our simulations (r < 1 AU) once the mass 
accretion from the remnant of the first core finishes because 
the massive remnant still exists and is accreating onto the 
disc, even at the end of the simulations. Unfortunately, it 
is almost impossible to investigate the further evolution of 
the disc without a sink. Although the sink may introduce 
numerical artefacts (especially in the few sink radius), it is 
an essential technique for investigating the long-term evo¬ 
lution of the disc. We will investigate the further evolution 
of the disc with the sink technique while remembering that 
this introduces numerical artefacts. 

In this paper, we showed that the SPH method 
is capable of treating MHD and non-ideal processes 
in realistic astrophysical simulations. Our results are 
largely consistent with those of the recent non-ideal 
RMHD simulations with the static -mesh-refinement code 
dTomida, Okuzumi fc MachidallioiSl ). Thus, our method is 



x(AU) 


Figure 10. Solid lines show the ratio of the sum of the centrifugal 
force and the pressure gradient force to the radial gravitational 

force, qi = | —— ^ ^ -1, as a function of the radius. Here, p and 

'll are the pressure and the gravitational potential, respectively. 
The dashed lines show the ratio of the centrifugal force to the 

radial gravitational force, q 2 = | ^ |. The red, green, and blue 
lines show the results of models 1, 2, and 3, respectively. The 
epochs are the same as in figure |8] 


reliable and can be used for astrophysical simulations. Be¬ 
cause the SPH method is relatively easily implemented and 
more flexible than static-mesh-refinement code, it can be 
used as an alternative method for many astrophysical prob¬ 
lems in which the magnetic field play the important role. 

In the simulations presented in this paper, several ap¬ 
proximations were adopted. The influences of these simpli¬ 
fications should be investigated in future studies. For ex¬ 
ample, we used a fixed dust grain size of a = 0.035 fim 
and a fixed the cosmic-ray ionization of ^cr = 10“^^ s“^. 
The latter is not good approximation for the dense region, 
p ~ 10“^^ g cm“°. We also used a simple rigidly rotating gas 
sphere as th e initial condition. As Joos. Hennebelle fc CiardU 
ll2012t) and iMachida. Inutsuka fc Matsumot3 ( 20141 ) have 
pointed out, the initial density profile and the magnetic field 
configuration strongly affect the size of the circumstellar 
discs. In future, we will investigate how the differences in 
the initial configuration affect the disc evolution. 
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Figure 11. Toomre’s Q value as a function of the radius in the 
X direction. The green and blue lines show the results of model 2 
and 3, respectively. The epochs are the same as in figurejS] 
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